#include <cusp/csr_matrix.h>
#include <cusp/dia_matrix.h>
#include <cusp/io/matrix_market.h>
#include <cusp/precond/spai.h>
#include <cusp/frobenius_norm.h>
#include <cusp/multiply.h>
#include <cusp/elementwise.h>
#include <cusp/detail/profiler.h>

#include <iostream>

int main(int argc, char** argv)
{
	typedef int IndexType;
	typedef float ValueType;

	cusp::csr_matrix<IndexType, ValueType, cusp::device_memory> A;
	cusp::csr_matrix<IndexType, ValueType, cusp::device_memory> P;
	cusp::dia_matrix<IndexType, ValueType, cusp::device_memory> I;

	if (argc < 2)
	{
		std::cout << "Missing argument." << std::endl;
		return 1;
	}
	std::cout << "Reading matrix from file: " << argv[1] << std::endl;
	cusp::io::read_matrix_market_file(A, std::string(argv[1]));

	std::cout << "Frobenius norm of A:    " << cusp::frobenius_norm(A) << std::endl;
	cusp::precond::spai<ValueType, cusp::device_memory>  M(A);

	cusp::multiply(A, M.matrix, P);
	cusp::gallery::eye(I, P.num_rows, P.num_cols);
	cusp::subtract(P, I, P);
	std::cout << "Frobenius norm of AM-I: " << cusp::frobenius_norm(P) << std::endl;
	cusp::io::write_matrix_market_file(M.matrix, "m.mtx");

	CUSP_PROFILE_DUMP();

	return 0;
}

